library(AlgDesign)
library(plyr)
library(Matrix)
#################################
# Step 1: Simulate data
#################################

ni <- 200 # num of individuals
np <- 5   # num of purchases for each individual
n <- ni * np
K <- 30  # num of brands
p <- 2   # number of individual specific vars

set.seed(11)
# data structure
sid <- rep(1:ni, rep(np, ni))
id <- rep(1:n, rep(K, n))
da <- data.frame(id = id, 
                 sid =  rep(sid, rep(K, n)),
                 choice = rep(c(1, rep(0, K - 1)), n),
                 brand = rep(1:K, n),                  
                 x1 = rnorm(n * K),
                 x2 = rnorm(n * K),
                 x3 = rnorm(n * K))   #brand specific vars
# individual specific variables
da2 <- matrix(rnorm(ni * p), ni, p)
da2 <- data.frame(sid = 1:ni, da2)
names(da2)[-1] <- paste0("x", 4:(3 + p))

#add in individual specific vars
da <- join(da, da2) 


# generate attributes
nG <- 5             #number of attributes
nlv <- 4            #number of levels in each attributes
Z <- gen.factorial(nlv, nG, varNames = paste("z", 1:nG, sep = ""))
Z <- Z[sample(1:nrow(Z), K), ]

da <- cbind(da, data.frame(z1 = factor(rep(Z[, 1], ni * np)),
                           z2 = factor(rep(Z[, 2], ni* np)),
                           z3 = factor(rep(Z[, 3], ni* np)),
                           z4 = factor(rep(Z[, 4], ni* np)),
                           z5 = factor(rep(Z[, 5], ni* np))))

da <- transform(da, id = factor(id), 
                brand = factor(brand), 
                sid = factor(sid))
cid <- rep(1:n, rep(K, n))

X <- model.matrix(~ z1 + z2 + z3 + z4 + z5 + x1 + x2, da)[, -1]
if (qr(X)$rank != ncol(X))
  stop("X does not have column rank")

# convert for object "mlogit.data"
da <- mlogit.data(da, choice = "choice", shape = "long", alt.var = "brand")
